Unified decoupling scheme for exchange and anisotropy contributions and 
temperature-dependent spectral properties of anisotropic spin systems 



R. Bastardis 1 , U. Atxitia 2 , 0. Chubykalo-Fesenko 2 , and H. Kachkachi 1 
1 PROMES-CNRS UPR8521, University de Perpignan via Domitia, 
Technosud, Rambla de la thermodynamique 66100 Perpignan, France and 
2 Instituto de Ciencia de Materiales de Madrid, 
CSIC, Cantoblanco, 28049 Madrid, Spain 
(Dated: March 6, 2013) 

We compute the temperature-dependent spin-wave spectrum and the magnetization for 
a spin system using the unified decoupling procedure for the high-order Green's functions 
for the exchange coupling and anisotropy, both in the classical and quantum case. Our 
approach allows us to establish a clear crossover between quantum-mechanical and classical 
methods by developing the classical analog of the quantum Green's function technique. The 
results are compared with the classical spectral density method and numerical modeling 
based on the stochastic Landau-Lifshitz equation and the Monte Carlo technique. As far as 
the critical temperature is concerned, there is a full agreement between the classical Green's 
functions technique and the classical spectral density method. However, the former method 
turns out to be more straightforward and more convenient than the latter because it avoids 
any a priori assumptions about the system's spectral density. The temperature-dependent 
exchange stiffness as a function of magnetization is investigated within different approaches. 

PACS numbers: 75.10.-b General theory and models of magnetic ordering - 75.30.Ds Spin waves - 
75.10.Jm Quantized spin models - 75.10.Hk Classical spin models 

I. INTRODUCTION 

Spin systems offer a rich area for fundamental research, always providing us with new open 
and challenging issues. In the context of modern applications, magnetic systems at the nanoscale 
have opened a huge laboratory for testing and applying the available methods with the challenge 
to adapt them to the constraints of the new area of magnetic nanoteclmology. Indeed, the rapid 
development of computers has opened a new trend for the magnetic materials design. Today the 
large scale materials modeling is often used as an efficient way to find optimal material performance 
in technological applications such as magnetic recording. The micromagnetic simulations represent 
now a powerful tool, especially after the development of publicly available software codes. To 
provide reliable predictions, the modeling methods should be improved with the incorporation of 
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detailed information from microscopic materials parameters into the macroscopic parameters such 
as magnetization, anisotropy or exchange stiffness. An additional problem arises when the full- 
fledged well-known approaches have to be extended to finite-size systems with acute boundary 
problems. 

Furthermore, multiple recent applications require temperature-dependent macroscopic proper- 
ties. These important ^applications include heat-assisted magnetic recording [lj], laser-induced mag- 



netization dynamics 
domain wall motion 
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2 1, thermally-assisted magnetic random memories [3fl and thermally-assisted 
4]. In this context the multi-scale scheme where the temperature-dependent 
macroscopic parameters are previously calculated numerically or analytically with the aim to use 
them in larger scale modeling has been proposed 

a a. 

The variety of methods, classical and 
quantum, analytical and numerical, were developed in the past and can be adjusted today for ap- 
plications within this multi-scale modeling framework. It is then necessary to take stock of the 
various methods, compare them and establish their respective limits of applicability. This is a 
tremendous task that has to be tackled before one can apply these methods to design new magnetic 
materials. 

Accordingly, the present work is about a few standard methods used for investigating the spec- 
trum of spin waves (SW) in magnetic systems at finite temperature and for arbitrary spin. These 
are the quantum Green's function (QGF) technique and its classical limit (CGF), the classical spec- 
tral density (CSD) method, and the purely numerical methods, i) one that consists in solving the 
stochastic Landau-Lifshitz equation (LLE) Q, \^ and ii) the Metropolis Monte Carlo (MC) method 
[lo| . Although these methods exist in the literature in different and multiple formulations, no sys- 
tematic comparison with the aim to establish their agreement and crossover has been made. One of 
our objectives here is to compare these methods and establish the best framework for the calculation 
of the temperature-dependent SW spectrum and physical observables such as the magnetization 
and susceptibility. For each method we discuss the most reliable implementation which gives the 
best agreement with numerical techniques and provide a clear crossover between the classical and 
the quantum case. In this task we have realized that no unified decoupling scheme used to take 
into account both the exchange and anisotropy contributions in the classical case has been given in 
the literature so far. However, this is exactly what is required for the purpose of the hierarchical 
multi-scale modeling, where classical Heisenberg-like Hamiltonian is parameterized via ab initio 
calculations and is used to evaluate temperature-dependent macroscopic properties 

aa 

. On the 

other hand, in the future the use of classical systems may be avoided if direct reliable calculations 
of macroscopic properties at the nanoscale based on quantum spin systems are available. This is 
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why it is important to establish clear connection between quantum and classical approaches. 

It is well known that the Green's function and spectral density methods involve a decoupling of 
high-order spin correlations into two-point correlations. Here we revisit this issue and demonstrate 
a clear connection between, on one hand, the classical and quantum approaches, and on the other, 
the CGF technique and the CSD method. In the quantum case, the spin operators satisfy the 
SU{2) Lie algebra and this implies that two spin operators commute when they refer to distinct 
lattice sites. In particular, the longitudinal and transverse spin fluctuations are uncorrelated when 
they refer to two distinct lattice sites and they are strongly correlated otherwise. However, a 
decoupling that may be successful in dealing with the exchange coupling contribution, at least 
at low temperature, may turn out to be a bad approximation for the local (with identical lattice 
sites) contributions of (on-site) anisotropy. This is why mean-field theory (MFT), random- phase 
approximation (RPA) and the Bogoliubov-Tyablikov approximation (BTA), which assume that the 
longitudinal and transverse fluctuations are uncorrelated, provide a reasonably good approximation 
for exchange whereas they provide rather poor results in the presence of anisotropy. 

Here we pay a special attention to this issue and considerably clarify the situation regarding 
the decoupling scheme that is used for exchange and anisotropy contributions. More precisely, we 
provide a unified decoupling scheme for both exchange and anisotropy contributions, for classical as 
well as quantum spins. Then, using this decoupling we obtain workable (semi-)analytical expressions 
for the SW dispersion, the magnetization and the critical temperature, which are supported by the 
good agreement with the numerical results of the LLE method and Monte Carlo (MC) simulations. 

In section |H] we define the generic system we study using the Dirac-Heisenberg Hamiltonian. In 
section HTT1 we discuss the various decoupling schemes used in QGF technique and show how they are 
related, and compute the SW dispersion and the magnetization. For the latter we also provide the 
analytical asymptotes at low temperature and near the critical point. Next, we work out the classical 
limit of this approach and obtain the corresponding dispersion and magnetization. For the latter, we 
pinpoint an interesting connection between the Callen's (quantum) expression for the magnetization 
and the MFT-like expression in terms of the Brillouin function in the quantum case. Apart from 
its elegance, this formulation makes it straightforward to derive the classical limit in terms of the 
Langevin function. We then turn to the CSD method and clarify the relevance of the decoupling 
scheme when it comes to treat the exchange coupling and anisotropy. We end this section with a 
brief account of the LLE and MC methods and a few expressions and numerical estimates of the 
critical temperature. Section IIV Al presents the results for the SW spectrum, the magnetization as 
a function temperature and field, and spin stiffness. The paper ends with a conclusion and outlook. 
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In Appendices [A] we present the main steps of the QGF for finite temperature, arbitrary spin, and 
oblique magnetic field. In appendix [Bj we give a detailed demonstration of some expressions used 
within the CSD approach. In appendix [0 we give a few expressions and numerical estimates for 
the critical temperature. 

II. MODEL: HAMILTONIAN AND SYSTEM STUDIED 

We study a spin system of M atomic spins Si = Ssi, with |sj| = 1, interacting via a nearest- 
neighbor exchange coupling Jy. In addition, each spin evolves in a (local) potential energy that 
comprises an on-site anisotropy and a Zeeman contribution. The anisotropy is taken uniaxial with a 
common easy axis pointing in the z direction; the magnetic field is applied in an arbitrary direction 

so that H = He^. The Hamiltonian of the system then reads 

M M 
^ = - 2 E J *i S * • S i " K H - (MbH) J> • e h . (1) 

(id) i=1 i=1 

We consider only box-shaped systems of size N x x N y x N z = N with, e.g. a simple cubic (sc) 
or a body-centered-cubic (bcc) lattice structure. 

III. METHODS 

We would like to investigate the spectral properties of such systems using and comparing two 
groups of methods: i) the (semi-)analytical methods, namely the classical spectral density method 
and the classical or quantum methods of Green's functions (GF) at finite temperature, ii) the 
numerical methods that consist either in solving the stochastic Landau-Lifshitz-Gilbert equation 
(LLE) in the Langevin approach or Monte Carlo (MC) simulations. 

A. Quantum Green's function approach 

The Green's function approach has been used thoroughly in almost all areas of physics. For spin 
systems, this approach allows us to obtain and investigate all kinds of observables. As compared 
to spin- wave theory (SWT), it makes it possible to obtain in a more systematic way the excitation 
spectrum at finite temperature for arbitrary atomic (nominal) spin. 

For our present purposes, we re-derive the basic equations involved in this approach and apply 
them to the Hamiltonian ([!]). In the latter the magnetic field is applied in an arbitrary direction 
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with respect to the (common) anisotropy easy axis and as such, a slight reformulation of the basic 
equations is needed with respect to the equilibrium configuration. In particular, Callen's formula 
11| for the magnetization in the case of arbitrary spin has to be re-derived in this context. The 
details of these calculations are given in Appendix 1X1 

We introduce the retarded many-body Green's functions 

gr{i-j,t) = g^(n - Tj ,t) = («(t);o7(0))) r = -iB(t) ([<(*), <(o)]> • (2) 

where <Tj are the new spin variables obtained after rotation of the original variables Sj to the system 
of coordinates where the z-axis coincides with the direction of the net magnetization [see Appendix 
IX] ; (...) denotes the usual thermal average. Then, one establishes the equations of motion of the 
GFs Gij~ , Gjj~ , Gfj~ whose solution renders the SW dispersion. 

The equation of motion for a GF of a given order in spin operators generates GFs of higher 
orders and this leads to an infinite hierarchy of GFs satisfying an open system of coupled equations. 
In order to close this system of equations and solve it (in Fourier space), one is led to apply a certain 
scheme for breaking high-order GFs into lower-order ones, thus adopting a certain approximation 
of the magnon-magnon interactions. Finding an adequate scheme for doing so has triggered many 
investigations each dealing with a specific situation with a particular Hamiltonian. Unfortunately, 
there is no general or systematic procedure. In fact, the variety of decoupling schemes only reflects 
the complexity of dealing with magnon-magnon interactions and thereby the nonlinear SW effects. 
In the following section, we present a discussion of the main decoupling schemes known in the 
literature and also propose some improvements that allow for a certain unification thereof. 

1. Decoupling schemes 

When applying mean-field theory (MFT), random-phase approximation (RPA), or the 
Bogoliubov-Tyablikov approximation (BTA), it is implicitly assumed that the longitudinal and 
transverse fluctuations are uncorrelated and this is a valid approximation only when they refer to 
distinct sites. Indeed, the idea behind this approximation consists in writing 



{[AB,C\)~{A)x([B,C\). (3) 

For spin systems, the factor (A) is usually the thermal average of <r 3 and thereby is related to the 
temperature-dependent magnetization. Hence, in practice one rearranges the various terms so that 
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<7 3 appears on the left and then use the approximation (|3|). However, in the (local) anisotropy 
contributions the product factors are at the same site and thus the longitudinal and transverse 
fluctuations are correlated, which turns this kind of decoupling schemes into relatively bad approx- 
imations. 

In Ref. [3| it was argued that one may avoid this approximation inherent to a decoupling 
scheme by establishing 25 equations of motion for the anisotropy functions. The problem with this 
approach, however, is that in practice one has to specify the spin S thus limiting the calculations 
to a particular material. In addition, it is not obvious how to obtain the classical limit from the 
final results. 

For the exchange coupling, RPA is commonly used with reasonable satisfaction since the cor- 
responding results for the SW dispersion and thereby the magnetization compare fairly well with 
other techniques such as Monte Carlo (MC) [see Ref. jl^ for a recent review], as long as the mag- 
netization curve at low temperature is concerned. However, for a more precise estimation of the 
critical temperature Tq, Callen's decoupling scheme turns out to be much more efficient, though it 
leads to a self-consistent equation for Tc which is more difficult to tackle analytically. Indeed, it 



was shown by Tahir-Kheli and Callen 



16l4l8| that the more sophisticated decoupling scheme 



af(r)a+(r);aj(0))) * (of) ((a^aj (0))) - (af) ^LZlI ((^(t);^" (0))) (4) 

takes, to some extent, account of magnon-magnon interactions and renders a nonlinear equation 
for the magnon dispersion w(k), see below. 

For on-site magneto-crystalline anisotropy the simplistic RPA decoupling leads to poor and 
even wrong results. In the presence of anisotropy with typical ratios K/J, the Anderson-Callen 
decoupling scheme, originally proposed by Anderson and Callen fisl . and later generalized by 
Schwieger et al. |l4l | to a rotated reference frame, turns out to be rather efficient in producing 
reasonable results. This is typically of the form 



K {{{afar + a'af) (t) ; aj (0))) « 2K a (af) x {{a~ (i) ; aj (0))) , (5) 
with the effective anisotropy factor 



IC a = K 



(6) 
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The identity 



fa^r) + <<W> = 2S i s + !) - 2 <^ 3 > . 



is derived from the quantum-mechanical identities 



<j~o + = 5(5 + 1)- (<t 3 ) 2 -ct 3 , 



a = - K a a 



(7) 



It is well known that the decoupling ([5]) is valid for all spin values 5 and renders good results 
when compared with the exact treatment of anisotropy and with quantum MC when K/ J is small 



Similar to the decoupling in Eq. ([5]), the following decoupling for anisotropy has been suggested 



KrH(r) + a!(r)at(r);aJ (0))) = « 
Then, splitting the right-hand side as follows 



25 2 
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^7(0) 
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((at(r);aj (0))) = <af> fl - ((^(r);^ (0) 



+ 
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we may propose the following decoupling 



((^(rK+(T);a7 (0))) = (af) 

Comparing this decoupling for the anisotropy contribution with Eq. ([4]) for the exchange con- 
tribution, we see that the former follows from the latter upon setting in the latter I = i, i.e., 
restricting the product of spin operators to the same lattice site. In fact, this is a consequence of 
the way af is written in powers of af and the products afaj. More precisely, if we start from the 
quantum-mechanical identities ([7]) and then multiply them by a and 1 — a respectively and add 
the resulting equations we obtain 

a 3 = aS(S + 1) - a (a 3 ) 2 + ~ 0^) a ~ a+ > ^ 



25 2 



(8) 



8 



where ais then determined so as to comply with the limits at zero temperature and near the critical 
point 



18|] . This leads to a = (a 3 ) /2S 2 . 
Next, we insert the expression ([9]) for cr 3 in products of spin operators such as those appearing 
on the left-hand side of Eqs. ((31 EJ) and use Wick's or RPA-like decoupling to obtain the decoupling 
for exchange and ([8]) for anisotropy contributions, respectively. In fact, there exist several other 
decoupling schemes in the literature with expressions for a that are polynomials of different degrees 
in m = (c 3 ) /S. Namely, a = corresponds to RPA (or BTA), a oc m to Callen's decoupling, 
a oc m 3 to the decoupling proposed by Copeland and Gersch (CG) an d 

aim) = lsWTT) m+ sWTT) m3 (10) 



to the decoupling proposed later by Swendsen 
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As already discussed, these polynomials with increasing degrees are approximations to the more 
rigorous calculation of spin correlations that consists in computing contributions of high-order of 
Feynman's spin diagrams as is done in Refs. 
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23l . As it will be seen later in section IIV Al the 



corresponding decoupling yields fairly precise results for the magnetization and critical temperature. 

2. Spin-wave dispersion 

Applying for instance the RPA decoupling to a homogeneous ferromagnet, i.e. with (erf) = (<7 3 ), 
(see details in Appendix [A]) one derives the magnon energy with respect to the equilibrium state 



£ 2 (k) = (Mk)) 2 = Al - Bl (11) 
= [{ gf x B ) (H x sin?? + H z cost?) + K. a (a 3 ) (2 cos 2 1? - sin 2 0) + J (a 3 ) (1 - 7k )] 2 
-(/C ff ) 2 (cj 3 > 2 sin 4 0. 

where [see Appendix [A] for notation] 

A = L 3 + /C (T <(7 3 >(2cos 2 0-sin 2 0)+Jo< ( 7 3 >(l-7 k ), 
#k = /C CT (cr 3 )sin 2 

J = J (0) being defined in Eq. (|A9|) . 

This dispersion relation is effectively obtained within the linear spin-wave theory, because the 
high-order GFs stemming from exchange contributions have been decoupled using the RPA which 
does not take account of spin correlations or magnon-magnon interactions. Indeed, following the 



standard procedure described in the appendix of Ref. Q (and references therein) one arrives at the 
equation for the dispersion (in the case $ = 0, i.e., of longitudinal field and (of) = {Sf}) 

hw (k) = fe k = fas) H + 2K S (S z ) + {S z ) [J(0) - J(k)] + £ [J(p) ~ J(P ~ k)] (np> , 

p 

(13) 

where JCs reads, 



K s = K a (0 = 0) = K 



1-^(S(S + 1)-(S Z S Z )) 



(14) 



and a = m for Callen's decoupling. (n p ) is the thermal occupation number given by the magnon 
Bose-Einstein distribution 

where /3 = 1/k^T. 

Then, using translation invariance we see that 

J(p) - J(p - k) _ 7 P - 7 P -k 



J(0) - J(k) 1 - 7k 



7 P 



and thereby 

fiwk = (sa*b) # + 2/C5 (5*) + Jo (S z ) Q (a, 0) (1 - 7k ) . (16) 
Here we have introduced the exchange "stiffness coefficient" 

Q^P) = 1 + is £ = 1 + 77s £ 7p <»p> ( 17 ) 

p p 

where a is, as defined earlier, depends on the decoupling scheme. 

3. Magnetization 

Now, we turn to compute the magnetization for an arbitrary spin S. In Callen's method [see 
Ref. [ll( and references therein] one considers the GF 



E^(i-j,t;0 = -i9(t) x ([af (t),exp(£aj(0)) <t?(0)]) 



(18) 
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Then, replacing the GFs in the system (|A15P by their analogs from Eq. (|18p we obtain the new 
system of EM 



-£ k E+ + (w + „4 k )~. 



<K + (0),exp(^f(0))a,r(0)]>^S+- 



0. 



Following again Callen's procedure we obtain for the first moment (er 3 



(19) 



<- 3 > 



{s -n)(i + n) 2S+1 + {s + i + fi)o 2S+1 



(i + nys+i - n2s+i 

where Q is the following function of the 1 st and 2 nd moments (c 3 ) , (cr 3 <7 3 )=C 33 , 



(20) 



0«a 3 ),C 33 ) 



1 1 

2 



-4k cQth / ffl^k 



1 



(21) 



A^ c is the number of unit cells in the Bravais lattice of the ferromagnet and is also the number of 
allowed wave-vectors in the Brillouin zone. 

Eq. (|20p . together with (|16p and (|2ip . constitutes a transcendental equation whose solution 
involves several sums (integrals) in Fourier space. In general, it is a heavy task to solve Eq. (l20l) 



especially for lattices with several sub-lattices. Nonetheless, it was shown in Ref. |24j | that the 
magnetization a = (er 3 ) in Eq. (f20|) can be recast into the following more compact form 



a = SB S [SX] 



where X is defined by 







- 1 

and Bs(x) is the Brillouin function (for the quantum spin S) 



(22) 



(23) 



B s (x) 



25 + 1 



coth 



25 + 1 
25 



1 / x 

25 COt H25 



(24) 



25 

More generally, it was shown that the higher moments {a 11 ) can all be expressed in terms 

of the reduced magnetization m = a/S whereby the temperature T and field H enter via m = 
m(T,H). It was argued that this model-independent MFT-like result stems from the exponential 
form of the probability density, i.e., p = e xs * /Tr e xsz . Indeed, Eq. (|23p expresses the fact that 
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in MFT all the excitations are degenerate and that one may define the energy e = X/k^T as the 
effective energy of the molecular-field-like excitations with the same occupation number as the true 
excitations. On the other hand, we note that Eq. (|22p is also a transcendental equation for a, 
similar to Eq. (|20p . though much more compact and it readily yields the classical limit, as will be 
seen below. In addition, this establishes the connection to the standard result of MFT. 

In order to solve either equation, i.e. (I20p or (|22|) . and obtain the magnetization m(T,H) one 
has to supplement the latter by a second equation for the correlation function C 33 = (V 3 <7 3 ); this 
is obtained by the Callen's procedure (for £ = 0, see Eq. (|18p ) which leads to 

{araf) = 2a x ft(a,C 33 ) . (25) 

together with 1 st the identity in Eq. ([7]). 
The latter also yields 

(a 3 a 3 ) = S{S + l)-a- (a~a + ) = S (S + 1) - (1 + 2fl) a. (26) 

Finally, the magnetization (in the rotated frame) a is given by the solution of the following 
system of two nonlinear (coupled) equations 

a = Si?s [Sin (1 + 1)], 

(27) 

C 33 = S{S + 1) - (l + 2«)ff. 

In general, this system can only be solved numerically as it involves transcendental equations 
with several integrals. However, we can establish a few analytical expressions for the magnetization 
in the limiting temperature regions T — > and T — > Tq and upon restricting ourselves to a 
longitudinal magnetic field, i.e. applied along the direction e3 (ip = 0). 

In this case the SW dispersion hu^ in Eq. (jlip simplifies into 

ftwk = Ak = L 3 + 2K, a a + J a (1 - 7k) 

and Eq. (|2ip becomes 

V ' J N c ^2 
c k 

Low temperature asymptote 



coth 



(28) 
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At low temperature, the spins are strongly correlated and thereby the correlation function C 33 
tends to S (S + 1). As a consequence, the effective anisotropy obtained from the Anderson-Callen 
decoupling scheme simply yields JCs — > K [see Eq. (|14p ] so that the system of equations (f2"T|) 
decouples leading to a closed equation for a whose solution then is Eq. (|20p . Expanding the latter 
in terms of il (which becomes small at low temperature), we find 

a ~ S - n (a) . 

Moreover, at low temperature only low-energy spin waves are excited and these are the long- wave 
length modes. Hence, in the limit of small wave vectors, we have the dispersion relation 

Huj k ~ g/_i B H z + 2Ka + Aak 2 , (29) 

where A = J5 2 , 5 2 = ^2 a 2 J (a) / ^2 J (a) , a is the lattice parameter and J(a) is the exchange 
coupling over the nearest-neighbor bond. 

Next, upon expanding f2 in terms of temperature T (or rather in k^T/ JqS) we obtain 



3r \ 



3/2 



'3/2 



h + kS 



(30) 



where 

oo 

Z p (x) = J2n- p e- nx , 
n=l 

and Z p (0) = £ (p) is the well known Riemann zeta function. We have also introduced the following 
dimensionless parameters 



t= J_ = ^bT h _ (g/i B ) H K _2K 
(3Jq Jo ' Jo Jo 

Obviously, in the present limit and in zero applied field, one obtains the well known "3/2" Bloch's 
power law for the thermal decrease of the magnetization. 

Near-critical temperature asymptote (H = 0) 

Just below the Curie temperature, in the absence of magnetic field, the mean number of excited 
quasi-particles and their density are large, and it is then a reasonable approximation to pass to the 
continuum limit. In this case, in Eq. ([28]) we make the transformation 



13 



where vq is the volume of the unit cell of the direct lattice. 

Next, in this limit the system (I27p again decouples and leads to the Callen's expression (I20p for 
the magnetization, similarly to the low-temperature limit. In addition, we may write for C 33 



C 



33 



5(5 + 1) 



(31) 



where the factor 1/3 stems from the three dimensional rotational symmetry (SO (3)) of spins that 
starts to recover as the temperature reaches the critical temperature of the ferromagnet. 

Consequently, upon inserting in Eq. (JT3D C 33 = (S Z S Z ) given by the result above, dropping the 
nonlinear SW contributions, and neglecting the second-order terms in a we obtain the dispersion 



hulk = 2Kr]a + a J (1 - 7k) = JqcfX 1 (1 - A7k) 



with 



77 = 1 



5(5 + 1) 



A = 



1 



(32) 



(33) 



35 2 ' '■ 1 + t/k" 

In addition, a is rather small because the density of SW is large and since Wk is proportional to 
a, as is seen in Eq. (f32|) . we can expand Q ((o -3 )) in powers of and obtain 



n (a) = v 



dk 



1 



(2vr) 3 eP^ - 1 V ° J (2vr) 3 \/3tiui k 
Let us now compute these integrals. Using ([32|) . the first contribution reads 



dk 



1 



1 (3?iwk 

2 12 



(34) 



dk 



AP(A) 



(2vr) 3 Phut 



-T. 



a 



where we have introduced the well known lattice Green's function [see Ref. 
therein] 



26 and references 



/ \ \ f dk 1 
Analytical expressions for this integra 



Ref. 



M, see also Eq. (4.2) in Ref. 



27 



(2vrri-A 7k - (35) 
for various limiting cases of the parameter A are given in 
In our case, from (|33l) we have A = (1 + Tjn)^ 1 ~ 1 — 7]K 



since k = 2K/J$ <C 1. Hence 5X = 1 — A <C 1 and according to Refs. 

P(A) ~ - c (1 - A) 1/2 . 



26 



27 we have 



(36) 
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W is the Watson integral that evaluates to 1.51639 for a sc lattice and to 1.39320 for a bcc 
lattice; cq is a lattice-dependent constant that is equal to - (f) 1 2 — 1.16955 for the sc lattice and 
to 2 3 / 2 /7r ~ 0.90032 for the bcc lattice 26]. Next, using the fact that for both sc and bcc lattices 
26| 

dk 



ill 



\2n+l 



(2vr) 3 

we compute the remaining contributions in (j34|) and obtain 



0, n = 0,1,2,... 



(37) 



AP(A)r 1 la 
a 2 + 12Ar' 



Finally, using this expression in Eq. (|20p and expanding with respect to a we obtain the following 
asymptote for the magnetization 



a 



QGF 



2V15AP (A) t 



3AP (A) r 



1 



(38) 



^/4S(S + 1) + 5P(A) - 3 V S(S + 1) 
This asymptotic expression is plotted in Fig. [2] where it favorably compares with the other 
(exact numerical) magnetization curves. 

Now, in this relatively high temperature regime, magnon-magnon interactions become relevant. 
In order to take them into account, we consider the dispersion in Eq. (|32p to which we add the last 
contribution in Eq. (fl6l) . i.e. 



tiwk = JoaQ(a, r)A 1 (1 - , 



where 



A 



K+Q(a,r) ■ 



(39) 



Q(a,r) = l + ^£ p -^. 



Upon neglecting the second-order terms in a the last expression leads to the transcendental 
equation for Q 



Q{a,r) 



1 + 



Q(a,r) 



W N -C Jl 



Q(<x,t) 
k + Q(a,r) 



(40) 



In the absence of anisotropy one can easily solve the latter and obtain 
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Q(a, r)~Ul + yi + 4^r(W0v-l)) = Qexch(a, r). (41) 

Since the anisotropy contribution is much smaller than that of exchange we may seek a solution 
for Q (a, r) in the form 



Q(a,r) ~ Q CX ch(a,r)(l + e), e = ^ exch . 

Vanis 

Then, inserting this in Eq. (|40p and expanding successively with respect to e and then with respect 
to k, we obtain (to first order) 



Q(ol,t) ~ <2cxch(a,r) - 4>— -J K . (42) 

2Qcxch(a, t) - 1 Y « + Qcxch(a, t) 

Next, using the same expansion for Q (^<r 3 )), similar to Eq. (I34D we get 

, AP(A) 1 la 

fi(cr) ~ 1 — 

v 7 cj 2 12Ar 

which leads to the following asymptote for the magnetization 



-QGF ^, v^5AP(A)r / n 3AP(A)r 



~ v ^ . / 1 v > (43) 

Q(a , T )J 5(5 + i ) + MAE3V S(S + l)Q(a, r) 



4 

Note that this expression reduces to that in Eq. (|38j) if we set a = since then Q(0,r) = 1 
and A = A, which corresponds to the RPA decoupling. On the other hand, as we will see later 
[see Eq. (IClj) ]. one has to use this expression instead of (|38]) to obtain the critical temperature. 
In addition, as far as the magnetization is concerned, Eq. (|43|) renders a more precise profile for 
relatively higher temperatures. 

B. Classical Green's function approach 

In many situations, the classical approach turns out to be appropriate for describing the magnetic 
properties of the system studied. Therefore, it is worth establishing analogous expressions as in 
the quantum case by carefully examining the corresponding decoupling schemes and controlling the 
various approximations. Accordingly, in this section we establish a complete procedure, analogous to 
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the quantum-mechanical one, that yields the classical SW dispersion and thereby the magnetization. 
In particular, we provide the classical analog of Callen's decoupling scheme, for both exchange and 
anisotropy contributions. 

For this purpose we first set $ = and return to the spin variables Sj. We then introduce 
the classical spin vectors s, = Si/S and make the substitutions Jy — > J'- = S 2 J{j,K — > K' = 
S 2 K, H — > H' = SH in the Hamiltonian ([I]). Next, we define the classical two-time (retarded) GF 

Ga(r) = ((si(r);sj(0))) = -<0(t)({«+(t), S T(0)}), (44) 
and its (time) Fourier transform 

oo 

Gijiu) = I drGij^e^ = ((*t u ■ ( 45 ) 

— oo 

Using the Poisson brackets for the classical spin variables Sj [3| 



i s fi s j} = ^ijsf, {«+ Sj } = -2iS ijS l 
we obtain the equation of motion for Gij(r) and thereby for its Fourier transform Gij(uj) 

- icoG tJ (u) = -2i5ij <sf(0)) - i {gn B H') G tJ {uj) - 2iK' ((sf (r)s+(r), S "(0)))^ 

+ *E J il ((«? W^W;*7(°))) - *E J * ((W4(t)]sJ(0))) . (46) 
l w i w 

Then, in analogy with the quantum-mechanical decoupling of exchange in Eq. @, we propose 

the following decoupling scheme 

(s + s~) 

((<(T).+(T);*r(0))) w * ^(f+W^^^-^^i^W^DL' ( 4 7) 

We will show below that this decoupling scheme leads to the correct classical limit of the SW 
dispersion and magnetization. 

Similar to Eq. ([8]), the decoupling of anisotropy contributions is obtained from the equation 
above upon setting I = i. Note that this way the same decoupling scheme applies to both quantum 
and classical spins, and to both exchange and anisotropy. As discussed earlier, for quantum spins 
this unification of exchange and anisotropy decoupling schemes is due to the expansion in Eq. ([9|) 



17 



for S z . However, on the classical side there is no such expansion. This is a consequence of the fact 
that the second identity in Eq. ([7]) becomes meaningless owing to [S + , S~] = 0. 

Therefore, applying these two decoupling schemes and passing to the Fourier space in Eq. (l46l) 
we obtain 



2M- 



m 



u>> - (J (k) 

with the classical dispersion relation (cj^ = ^ (k)) 



Me = 9VbH' + 2K'm 



+ m [ J'(0) - J'(k)] 



+ W^[ J ' (p) - J,(p - k) ]^ eiP ' riJ (^ T 



Note that we have used the translational invariance to write (sf) = (sf) = m. 



Now if we apply the classical analog of the spectral theorem 



29 



30(| , i.e., 



Gk{oj' + ie) — Gk{oj' — ie) = —AmNm5(uj' — uj' k ), 

we obtain 



Inserting these expressions back into we obtain the classical analog of the dispersion relation 
that accounts for the SW interactions 



huj f k = g^sH 1 + 2K'm 

+m [J'(0) - J'(k)] + 



rri 

2 

m 

w 



E — 



p p 



(49) 



E 



j/(p)_ J'(p-k) 
M 



We stress again that only after solving this transcendental equation, one obtains the final SW 
dispersion lj^. This is, however, a heavy procedure because also enters the magnetization m, 
which in turn involves via ^, an d vice versa. At each step one has to compute three-dimensional 
sums (or integrals) in Fourier space. 
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Obviously, this dispersion can also be obtained by taking the classical limit of the quantum GF 
result, i.e. Eq. (jl3[> . Indeed, in the presence of uniaxial anisotropy, the Anderson-Callen decoupling 
yields the equation [see Ref. [3] and references therein] 

1 



= g^ B H + 2K(S z 



-^(S(S + 1)-(S Z S Z )) 



+ (S*) ( J(0) - J(k)) + £ 



J(p) - J(p - k) 

e (5hw p _ i 



Then, using the identities ([7]) and making the substitutions uj^ = u'^/S, Jij = Jlj/S 2 ,K 
K'/S 2 , H = H'/S, together with m = (S z ) fS, we obtain 



Huj' k = gn B H' + 2K' 



m 



m 

l + 2S 



{8+8-) 



+ m(j'(0)-J'(k))+^£ 



j/(p) _ j/( p _ k ) 
- 1 



In the classical limit (n p ) in Eq. (|15p becomes 

1 



< n P> 



hull 
g IE 



(50) 



and thereby 



S^oo 



gH B H' + 2K'm 



(s+8-) 



+ m[j'(0)-J'(k)]+^£ 



■/'(p)-J'(p-k) 



Next, upon replacing (s^Sj ) by its expression given in Eq. (148|) we obtain 

+ m [/(o)-/(k)]+^x; 



m x - 1 

^ p V P/5->oo 



(51) 



J'(p)-J'(p-k) 



5->oo 



This is the dispersion in Eq. (I49p . which was obtained directly from the retarded classical GF (1441) 
using the (classical) decoupling scheme (j4"7|) for exchange and its analog for anisotropy. Therefore, 
starting directly with GFs for classical spins and using the classical analog of the spectral theorem 
leads, as it should, to the same result that is achieved by proceeding with the GFs for quantum 
spins and taking the classical limit at the very end. 

Similarly to the quantum case, the dispersion (|5ip can be recast in the form 



S— >oo 



gii B H' + 2/C'm + mJ'(0)Q' (a, f3) (1 - 7k ) 
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where we have introduced the classical analogs of the effective anisotropy (|14p and the exchange 
stiffness ([Tj 



K! = K' 



/37V ^ 



1 



7 P 



p ^ \ VIS— >oo 



Before ending this section we discuss the magnetization. The large-spin limit, i.e. S 
yields the classical limit of the Brillouin function, that is the Langevin function, i.e., 



oo, 



lim Bs (x) = C (x) = coth (x) . 

S^-oo x 
On the other hand, this is what one obtains when the quantum spins are replaced by classical 
vectors and, in the partition function, the trace operator is replaced by integrals on the spin variables 
(or their spherical coordinates). Doing so for independent spins in a magnetic field x leads to the 
Langevin function. 

Now, in Eq. (122|) setting m = (c 3 ) /S and taking the limit 5 — > oo yields the magnetization in 
the classical limit, i.e. 

lim m = < S *) class =C(\), (52) 

S—s-oo VP/ 

with 



being the classical density of SW excitations. 

We note in passing that it is more straightforward to obtain the classical limit (|52[) from Eq. 
(|22p than from (|20p . On the other hand, Eq. (|22[) provides a clear connection with MFT. Indeed, 
as discussed earlier, this connection can be revealed by noting that all quasi-particle excitations in 
MFT are degenerate and thus one can simply drop the dependence on the wave vector in Eq. (|28D . 
However, this similarity in form should not shadow the fundamental difference, namely that in pure 
MFT the magnetization (<s 2 ) c i ass is calculated self-consistently using (in a longitudinal magnetic 
field) 

(s z ) = C [f3S {gfi B H' z + 2K' (s z ) + J' (0) (s z ))] (54) 

while in Eq. (|52l) one explicitly takes into account the SW dispersion via p. This SW den- 
sity is obtained by the GF technique using the RPA decoupling for exchange contribution and 



20 



the Anderson- C alien decoupling for single-ion anisotropy contribution. Eq. (|52|) is also a (self- 
consistent) transcendental equation because p is a function of the dispersion u^. For a comparison 
of the corresponding critical temperatures see Appendix O 

Magnetization asymptotes 

The low-temperature asymptote for the magnetization is obtained by expanding p and then the 
magnetization with respect to r = k-gT '/ Jq. Neglecting the terms due to Callen's decoupling for 
exchange and anisotropy terms in the dispersion relation (l49l) we obtain 

ftj^ = gp B H' + 2K'm + m ( J'(0) - J'(k)) . 
which may be rewritten as (J'(0) = Jq = zJ') 



huj' k = Jq [h' + (1 + k) to] (1 - i/j (to) j k ) 
where we have introduced the function 

ip (m) 



with 



ti 



h' + (1 + k) m 



9HbH' _ gpBH _ h 
~ S 



(55) 



4 



JoS 



(56) 



and 2K'/J = 2K/J = k. 

Then, the density p in Eq. (|53p becomes 

T ' 

p = ip (m) Pn W x — 

m 

where r 1 = t/S 2 and the function Pjy [ip (m)] reads 

which is the analog of (I35p for a finite lattice of linear size iV. Asymptotic expressions of the lattice 
Green function P/v (G) without the zero mode (k = 0), for free boundary conditions (fbc) and 



31|, 



32j. Wjv is the well known lattice 



periodic boundary conditions (pbc), can be found in Refs. 
sum whose large-size (continuous) limit is the Watson integral W, introduced earlier in Eq. (j36|) . 
For a bcc lattice we have 

W N ~ W hcc (1 - , for fbc 
Wiv ^ W bcc (1 - 5j83) , for pbc 



21 



and for a sc lattice [32fl 



W N * W sc + 9 -*^, forfbc 



^ W N ~ W*. (1 - 5j9&) for pbc. 
W sc and Wbcc are the Watson integrals for the corresponding lattices and are given after Eq. 



Now, at low temperature the density of SW p is small and using C (x) ~ 1 — 1/x for large x in 
Eq. (I52p we obtain the asymptote for the magnetization (up to 2 nd order in r), upon expanding p 
around m — 1, 



2 / /\2 



(57) 



Note that ^ (1) is a function of the applied field, since according to Eq. 055 p . ip (1) = 1/ (/i' + 1 + k). 
To first order in r' Eq. (i57l) obviously recovers the low-temperature linear decay of the magnetization 
which is typical of the classical Dirac-Heisenberg models. At very low temperature we can neglect 
the second-order terms and expand with respect to the field h' leading to 



CGF 1 w n , . c h' + K , 

m ~ 1 r H \ r . 

1 + k 1 + k\ 1 + K 

This is also the SW theory result obtained in Ref. 



32j, second line of Eq. (65), in the absence 



of anisotropy. We remark in passing that in this reference SW theory was extended to account for 
finite-size effects in fine magnetic particles. One of the consequences of these effects is that there 
appears a critical field Hy ~ T/M, that corresponds to the suppression of the global rotation of 
the particle's net magnetic moment, below which the magnetization is quadratic in the applied 
magnetic field. In the present work, the system size is big enough so that Hy vanishes and the 
quadratic behavior of the magnetization is suppressed. A more thorough comparison of the present 
work with that of Ref. [32j will be addressed in a future work. 



Near the critical temperature and in the absence of the applied field we have 

p = \'P (A') -. (58) 



We note that we have replaced the finite-sum lattice Green function Pjv (^') by its continuum limit 
defined in (f35|) as this is appropriate in the present high temperature regime. In the absence of the 
magnetic field and neglecting Callen's decoupling for anisotropy and exchange, i.e. for h! = 0, we 
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have ip(m) = 1/ (1 + k). Likewise, i] in Eq. (|33p is simply replaced by one and thereby ip equals 
the parameter A but here with the "primed" parameters, i.e. A — >■ A' = 1/ (1 + k). 

Then, since the magnetization m is small the density of SW p is large. Hence, using C (x) ~ 
x/3 — x 3 /45, for small x, and solving for m CGF we obtain the asymptotic expression 



m 



CGF 



V^A'P (A') r' \J\ - 3X'P (A') t'. (59) 



As in the quantum case, it is possible to take into account the magnon-magnon interactions 
from the last term in Eq. (|49j) . The magnetization is then given by the following expression 



CGF ^AP(A)rW 3AP(A)^ 

with a = 0, m 2 , m 4 , m (m + m 3 ) for the RPA, Callen, Copeland-Gersch or Swendsen decoupling, 
respectively. Now we can see that this can be recovered, as it should, as the classical limit of 
the asymptote (|43p obtained for quantum spins. Indeed replacing the various parameters by their 
classical counterparts {e.g. Jq by J' Q = JqS 2 , t' = t/S 2 , etc.) and dividing a by S we obtain 



^QGF 7l5AP(A)r' / 3AP(A)r' 



This readily yields the asymptote in Eq. (|59p upon taking the limit S — > oo, here and in Eq. 
39|), and writing m CGF = limg^oo (a- QGF /5). 



331 ] where 



We note that while submitting the present work we became aware of a recent work 
the classical GF method is developed along the procedure employed by Callen for quantum spins 
based on the generalized GF in Eq. (fl~8|) . The results obtained by the authors for the dispersion 
and magnetization are quite similar to ours. We stress, however, that the classical GF method we 
develop here is more straightforward as it avoids the difficult algebra involved in the calculation of 
the GF (fT8|) . which was introduced for dealing with arbitrary quantum spin S fll| . Moreover, our 
approach is based on a unified decoupling scheme for both exchange and anisotropy contributions 
and establishes a clear connection with the quantum-mechanical Callen's decoupling. In fact, the 
work in Ref. [33] about Callen's method together with the present approach provide a complete 



picture of the GF technique for classical spins. 
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C. Classical spectral density method 

In this section we summarize the basic ideas and formulas of the classical analog of the spectral- 
density method, the so-called classical spectral density method (CSD). One of the objectives of this 
method is to provide systematic and non trivial approximations in classical statistical physics when 
applied to classical spin systems. To the best of our knowledge, this was initially formulated in 3^| 
and later developed and applied by several authors [see Refs. Q,Q and references therein]. This 
approach is then compared to the classical GF technique developed in the previous section. In Ref. 
[33j|, the CSD method was compared to the classical analog of Callen's method. 

Here the spin Sj is a classical vector and the magnetization is defined by fh = (s z ). One 
then defines the classical spectral density (cj) of the time-dependent spin correlations. Then, the 
calculations proceed by assuming a given form (e.g. a Gaussian or a Lorentzian) for (w) involving 
some parameters The latter are obtained by solving a hierarchy of (moment) equations which 
are in turn obtained from a chain of equations for Green's functions of all orders f^ . For the 
Hamiltonian H in Eq. ([1]) one obtains the following dispersion relation (3, Q] 

Mc = h' + K- 2k ' + 4 - J k-q) (4 S = q> + 2 i 2k> + J Q " J k-q) <«^Iq>] ■ (61) 



q 



This involves the two correlation functions (si- s -k) anc ^ ( s £l-k) which have to be dealt with in 
order to proceed any further. (sj^s~ k )> is easily obtained as [7|| 

The second approximation made in CSD - the first being the form chosen for the spectral 
density A^ (u))- concerns the unavoidable decoupling scheme that is required for the calculation of 
the longitudinal correlation function (s^-k) - Here we stress that, as is seen from Eq. (f6"T|) . this 
contribution stems from the exchange as well as the anisotropy contribution. However, as discussed 
earlier, the decoupling that should be applied to the one or to the other contribution is rather 
different from the physical point of view since this depends on whether this contribution is a local 
or a bi-local term. Hence, let us summarize the results of our developments concerning this issue, 
which is extremely important as the soundness of the results is strongly dependent on its outcome. 

In fact, to obtain a decoupling for (s k sf_ k ) that stems from the exchange contribution, we 
may start from Eq. (|46p use the decoupling (|47p . and then Fourier transform the result. These 
developments are carried out in Appendix [B] and their outcome is the following exchange decoupling 
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scheme 



i 



\ S q/ \ S -q/ 



1 - m 2 ) <s+s_ q > 



(62) 



q q 

These calculations provide a clear derivation of the exchange decoupling used in the literature, see 



e.g. Ref. 



3d . Note, however, that the decoupling of the longitudinal correlation function (|62p is 



only valid under the sum over the wave vector q. 

For the anisotropy contribution one may start from the decoupling scheme proposed in Eq. (j47 
with j = i and decouple the high-order contributions, i.e., 

(a+8~) 



((^ k M;Vk(°))) [ 



i 



-)((-i k W;-q + k) 



In Ref. 



351 ] the higher-order spectral density was reduced to a simple form that leads to the 



correct results in the low- and high-temperature limits. This yields the following equation 



m / 1 



(8+8-) 



1 - - (s + ,<T 
2 X 



(63) 



In turn this renders the following expression for the magnetization 



with p being the spectral density defined in Eq. fj53|> . 

Here a remark is in order concerning CSD as compared to CGF. For a longitudinal magnetic 
field, Callen's expression (|20p or (|22p for the magnetization is exact, whereas expression (|64|) ren- 
dered by CSD is an approximation. Hence, in addition to the common approximation related 
with the decoupling scheme and which yields the SW dispersion, CSD introduces an additional 
approximation for the magnetization itself. 

The classical analog of ([7]) is obtained by using the condition |s| = 1 and the fact that for 
classical spins we have s + s~ = s~s + . That is 

( a V) = \-(s + 8-). 

Consequently, Eq. (|63p can be rewritten as 



2 A (s + s-)\ 

m v — —) 



= (s z s z )- l -(s + s~) 
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leading to the following decoupling for the anisotropy contribution 

(4 S -q) « m " + \ i 1 ~ m ") « S =q) ■ ( 65 ) 

As stressed earlier, we see that for the same longitudinal correlation (sqsf_ q ) we have a different 
decoupling scheme according to whether this results from exchange or anisotropy. Notice the 
difference in sign between Eq. (|62p and Eq. (|65p . 

Applying the decoupling (|62[) for the exchange and f)65[) for the anisotropy contributions to Eq. 
(|6ip we obtain the expression for the dispersion that coincides with the classical limit in Eq. ([49D . 
Summarizing, we see that only upon clearly identifying the origin (exchange or anisotropy) of the 
correlation function and applying the right decoupling scheme does one show that the CSD method 
renders the same results as the CGF technique. Next, we deal with the magnetization. 

Low temperature asymptote 

Expanding Eq. (fM|) with respect to p which is small here, we obtain 

-i ^ 2 2 

m ~ 1 — pm — —p m 

and then using the expression (|56p for p we get 

m ~ 1 - ^ (1) Pn [V> (1)] t - | (1) P N bp (I)]) 2 r 2 . 

Upon setting m ~ 1 in the right-hand side we see that this expression and the corresponding 
CGF asymptote (|5"7|) differ only at the second order in r' by a factor of 3/2. In the case ft' = we 
obtain 



m csD ^ j _ X , P (V) T ' _ | ( A 'p (V) T ') 2 . 
Near-critical temperature asymptote (/?/ = 0) 

Starting again from Eq. (f6~i|) in the absence of magnetic field and using the expression ([581 ) for 
p, we obtain 



m CSD ^ 



(66) 
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Similarly to the CGF method, if we take into account the magnon-magnon interactions by 
introducing the parameter a, the magnetization becomes 



m CSD ~ / / 1 

For comparison, we give the following relation between the CGF and CSD high-temperature 
asymptotes 



1 / 3AP (A) t 1 



CGF ^ yi5AP(A)T> AP(A)r' CSD 
" Q'(«,r') V ^-') ' 

D. Numerical methods 

One of the numerical techniques used here is based on the Langevin dynamics simulations of 
thermally excited spin waves [(], Q] in the classical case. The method is based on the numerical 
integration of the stochastic LLE 

| = -J(s t xHf + A S! x x Hf ] ) (67) 

where s is the classical localized spin corresponding to a localized magnetic moment with modulus 
fi. A and 7 are the Gilbert damping parameter and the gyromagnetic ratio respectively. The 
effective field, Hf ff , is then given by: 

Hf^W-if. (68) 

Here i s the stochastic term that describes the coupling to the external heat bath. Thermal 

fluctuations are included as a white noise term (uncorrelated in time) which is added into the 
effective field. The thermal fields are calculated by generating Gaussian random numbers and 
multiplying by the strength of the noise process. The correlators of different components of this 
field are given by 



(Ci,«(*)Ci,0(O> = ^^<%<M(* - ( 69 ) 

where a, {3 refer to the Cartesian components and T is the temperature of the heat bath to which 
the spin is coupled. 
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Using this technique, we simulate a generic three dimensional ferromagnet with a Heisenberg 
Hamiltonian as in Eq. ([1]) with an external applied held H parallel to the z— axis. The correlated 
magnetization fluctuations introduced by the random Langevin held are dealt with by Fourier 
analysis, both in space and time. More precisely, we transform the magnetization fluctuations 
s (r, t) = (s x (r, t) , s y (r, t)) around the equilibrium direction along the axis z via a Discrete Fourier 
Transform T>TT ', 

s{k,u n )=VFT(s(r,t n )) (70) 
where {t n } are discrete time instants and the wave vector for a hnite box-shaped ferromagnet with 



periodic boundary conditions takes the form 



31 



36l | ak a = 2Trn a /N a with n a = 0, 1, . . . , N a — 1; 
a = x,y,z. Then we compute the power spectrum density defined by -F(k, w) = |s(k, w)| 2 . 

The second numerical method used in this work is the classical Monte Carlo simulation technique 
using the standard Metropolis algorithm, see e.g. Refs. 



10. 



32j. The results of this method are 



used as a benchmark for those rendered by the (semi-) analytical methods of QGF/CGF and CSD 
with various decoupling schemes. For equilibrium properties it is well known that MC and LLE 
render similar results, with the difference that the former method is computationally faster at high 
temperatures whereas at low temperature LLE is faster. At the same time, we should note that 
the MC techniques do not include proper magnetization dynamics and thus are not suitable for the 
calculations of the spin wave spectrum but certainly are for the magnetization. 



IV. RESULTS AND FURTHER COMPARISON BETWEEN DIFFERENT METHODS. 

In this section we present a sample of the results for the SW spectrum and magnetization as a 
function of temperature and magnetic field, taking account of magnon-magnon interactions within 
various decoupling schemes. The second objective here is to compare the latter and assess their 
validity. We also evaluate the temperature-dependent exchange stiffness and provide (in Appendix 
ICj) analytical expressions for the Curie temperature within the decoupling schemes considered. 



A. Temperature-dependent magnetization within different decoupling schemes 

First, as an illustration of the temperature dependence of the SW spectrum, we plot in Fig. [T] 
the dispersion as a function of the wave vector k along the z axis, for different temperatures. 
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Figure 1: Dispersion relations, obtained by the CGF method, with wave vector (0,0, k) and without magnetic 



It can be seen that u^, which includes magnon-magnon interactions, is strongly dependent on 
temperature. At temperatures near the critical value, the SW softening is clearly seen. A favorable 
comparison of these curves obtained by the CSD method with those rendered by the numerical LLE 
method was presented in Ref. Q. 

Now, we present the magnetization curves, as a function of temperature and applied field, 
computed with the different methods for the bec lattice and iron parameters (per atom) J = 
1.44 1(T 21 J and K = 5.4 1(T 24 J. 

In Fig. [2j we plot the magnetization m = (S z ) /S as a function of (reduced) temperature 
r = k^T/Jo in zero magnetic field, as obtained from i) QGF with two values of the nominal spin 
S = 5/2,30, from ii) CGF, and from iii) classical MFT (CMFT), i.e. Eq. ((54)) ■ We see that as S 
increases the magnetization curve m tends to that rendered by CGF and CMFT. In particular, at 
low temperature we do see the evolution from the m ~ T 3 / 2 Bloch law to the linear law m ~ T, as is 
typical of the classical Dirac-Heisenberg model. It is interesting how the CMFT result agrees with 
that of CGS when m is plotted against the reduced temperature. The low-temperature asymptote 
(|30p shows a good agreement with the QGF curve for T < Tq/4. Similarly, the asymptote in the 
critical region (i38l) also reproduces correctly the QGF curve. 



field. 
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Figure 2: Magnetization against reduced temperature. Comparison of i) QGF for S = 5/2, 30, ii) CGF, and 
iii) classical MFT. We also show the low-temperature asymptote (|3"0|) and the asymptote near the critical 
temperature (|38[) . The critical temperature r c is that of the method used for obtaining the corresponding 
magnetization curve and r = ksT / Jq. 

In Fig. [3l we compare the magnetization curves rendered by CGF [see Eq. (|52p ] for various 
decoupling schemes, with MC as a benchmark. Here, we prefer to plot the magnetization against 
the absolute temperature r so as to see how different are the critical values of temperature rendered 
by the various decoupling methods. 

It is seen that the decoupling schemes of Callen and Swendsen agree quite well with MC. On the 
other hand, Copeland & Gersch (CG) and RPA decoupling schemes render nearly the same curve 
m (r) that goes below the previous curves at high temperature. This is simply due to the fact that 
decoupling schemes with terms of high powers of m, e.g. 3 in the CG decoupling and in the second 
term in Swendsen's decoupling [see Eq. QIC)]) ], lead to a negligible contribution at temperatures 
nearing the critical value. On the contrary, contributions that are linear in m in the decoupling 
schemes, such as Callen's and Swendsen's, do improve the magnetization curve at all temperatures. 
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Figure 3: Magnetization curves rendered by different decoupling schemes. The methods are compared to 
Monte Carlo. 



In Fig. 31 we compare the magnetization rendered by i) CGF and its Langevin function in Eq. 
(|52j) and ii) CSD given by Eq. (I64p . within RPA, and the two results are compared to MC. Globally, 
CGF renders a magnetization curve that keeps closer to MC than CSD method, which does so only 
at low temperature and near the critical temperature, 
m 
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Figure 4: Magnetization obtained by CGF and CSD within RPA. The two methods are compared to Monte 
Carlo. Inset: difference between CGF, CSD (with RPA decoupling) and MC. 



In the inset we plot the three differences between the CSD, CGF, and MC. It is seen that large 
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deviations occur for T/Tq > 0.4. 

In Fig. 0we compare, for the simpler case of an sc lattice, the classical Green's function method 
with three decoupling schemes, with the numerical LLE method. It is seen that LLE compares quite 
well with CGF using the Swendson decoupling in almost the whole range of temperature. Note 
however, that in the numerical LLE method the finite-size effects are clearly seen in the critical 
temperature region, as is also the case with MC, while the analytical methods do not ignore such 
effects for they implicitly consider an infinite lattice. 




Figure 5: Comparison between the GF methods, the CSD method with RPA decoupling, the LLE approach 
and the classical Green's function for an SC lattice structure. 



Next, we study the relative magnetization variation 



~ l(Tifl) = m^0-mM 
V ' ' l-m(T,0) 

as a function of the applied field for different values of temperature, without anisotropy. The results 
are shown in Fig. [6j 

In the quantum-mechanical case, we may use the low-temperature asymptote f|3U|) to get (in the 
absence of anisotropy) 



6m (T, H) ~ 1 - -7^3/2 ( 9J ^] ■ 
y ' C (3/2) 3/2 V kBT J 

This can also be seen within the quantum linear SW theory which renders exactly the same 
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expression. It is clear from the behavior of Z v (x) that in the low-temperature regime 8m decreases 
when the temperature increases. 

In the classical case and at low temperature, we use the asymptote (|57p and obtain 

Sm CGF ~ 1 - (1 + k) + (1 + k) (A' - V(l)) P N (1> (1)) x r'. 

Here we see that this expression increases when the temperature increases since A' = 1/(1 + k)> 
if) (1) = l/(h' + (l + «)) for any ti > 0. 

However, it remains unclear why in the quantum-mechanical case there is a change of behavior 
at a particular temperature because it is difficult to derive an (approximate) analytical expression 
for the latter. Indeed, this would at least require to derive the magnetization asymptote in the 
critical region in finite magnetic field which, unfortunately, leads to a rather cumbersome expression. 
Nevertheless, quantum spin effects are attenuated at high temperatures and as such the quantum 
approach renders the same behavior for the magnetization as the classical one. 




Figure 6: Field dependence of the relative magnetization variation 8m, (left) from QGF with RPA decoupling 
and (right) from CGF (Callen (a) and RPA (b) decoupling) and MC. 

B. SW spectrum and exchange stiffness at finite temperature 

Now we discuss the exchange stiffness as a function of the magnetization taking account of 
nonlinear SW effects. As we have seen, taking account of these effects (or magnon-magnon in- 
teractions), through the various decoupling schemes, leads to a temperature-dependent dispersion. 
This dependence on temperature comes about through the magnetization m. Let us consider, for 
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simplicity, the case with the sole contribution from exchange coupling. Hence, we may define the 
SW stiffness D as follows 

hu' k = D (m) (1 - 7k ) , (71) 

assuming that all SW nonlinear effects [see the last term in e.g. Eqs. (TT3l I49j) ] arc booked into the 
function D (m). On the other hand, in the absence of applied field and anisotropy, and for a given 
decoupling scheme with the parameter a introduced in Eq. ([9]) we deduce from Eq. (|17p that 

D (m) = JqItiQ' (a, m) 

with 



P 

in the classical limit. 

In the general discussed in the previous sections, the dispersion to>k and the magnetization 

m are solved for by using the system of coupled equations and then D (m) is obtained by fitting 
the curves of u; as a function of the wave vector k in a given direction in Fourier space. Next, we 
substitute Hlj'^ = J^mQ' (1 — 7k) in Q' to obtain 

„/ / 1 aim v 7r> , a (W — 1) . , , 

Q'(a,/3) N ^l-7k' mQ'(a,/3) 

where W is the Watson integral for the given lattice. 

At low temperature, the magnetization is given by (CGF or CSD) 



1 W ' 

Q'(a,(3) T 

leading to 

''-(l-m) 
Then, when this is substituted in Eq. (|72p yields 

Q =1 + — 777- x — 

W m 

and thereby the spin stiffness D (m) reads 

D (m) = Jq?ti 



W-l (l-m)a 

1 H — X V 



(73) 



34 



Now, defining (f> = a (m) /m and q = (W — 1) /W, we write 

D (m) = JqIu [1 + <;<f> (m) x (1 — m)] . 
Finally, considering the fact that at low temperature 1 — m is small so that we may write 

D (m) ~ 4m [1 - (1 - m)]~^ (m) = m 1- ^^. 

a (m) and thereby (m) is given according to the RPA, Callen's, Copeland and Gersch or 
Swendsen decoupling scheme, see Eq. fjlOjl ei seg. For the RPA decoupling, for instance, a = 
and thus D (m) ~ m, as it should. For Callen's decoupling, a (m) = m leading to D (m) ~ m 1_< \ 
For a decoupling scheme with a (m) = am + 6?n 3 we make an expansion around m — 1 and obtain 
I? (m) ~ m [1 — ^ (a + b) (m — 1)]. 

In Fig. [7] we plot the exchange stiffness as obtained numerically from Eq. (|7ip and Eqs. (fT3l 
l49l) . for the decoupling schemes discussed in Fig. [3l 




Figure 7: Exchange stiffness against the magnetization for different decoupling schemes obtained by CGF. 

As is seen in Eq. (|73l) and confirmed by the numerical results in Fig. [7j the exchange stiffness 
depends on the decoupling scheme or the way the spin-spin correlations are tackled, especially at 
moderate temperatures. Obviously, the curves corresponding to the four decoupling schemes merge 
for m ~ 1 (very low temperature) and m ~ 0, i.e. at high temperature where they exhibit a linear 
behavior. 
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V. CONCLUSION 

We have established clear connections between the quantum/classical Green function technique 
and the classical spectral density method, and have compared them with the numerical methods 
of Monte Carlo and Landau-Lifshitz-Langevin dynamics. We have proposed a unified decoupling 
scheme for both anisotropy and exchange contributions, for classical as well as for quantum spins 
which allows us to establish a clear connection between the various methods and to obtain reasonable 
results for the magnetization and critical temperature. We have computed the spin-wave spectrum 
at finite temperature and inferred the magnetization as a function of temperature and field and 
have obtained the exchange stiffness, for various decoupling schemes. Asymptotic expressions for 
the magnetization have been given at low temperature and in the critical region, both for classical 
and quantum spins, and the crossover between them has been established. As far as the (semi- 
analytical methods are concerned, it turns out that the classical Green's function technique is 
more straightforward and does not require any a priori assumptions about the the system's spectral 
density. In particular, Callen's famous formula for the magnetization is recast in a compact form 
using the Brillouin function. This makes it straightforward to obtain the classical limit leading to the 
familiar Langevin function for the magnetization. However, the outcome is still a transcendental 
equation involving the spin-wave density, unlike the Langevin form one obtains from mean-field 
theory. 

In future work, we would like to extend the present calculations and the Green's function tech- 
nique to finite-size systems by taking account of boundary and surface effects, similarly to what 



has been done in Refs. 
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32l |. This should be useful for studying the dynamics of multi-layered 



magnetic systems and magnetic nanostructures. 
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Appendix A: Quantum Green Function Method 



In this section we briefly describe the quantum Green function method (QGM) in the case of 
an oblique magnetic field. 



3G 



In order to compute the spin-wave (SW) spectrum and the magnetization, one deals with the spin 
fluctuations with respect to the equilibrium configuration, which has to be determined beforehand. 
In practice, one assumes that there exists a net direction of the system's magnetization denoted by 
e 3 

m = — r y^Si = me 3 . 

i 

We start by passing to the new coordinate system in which the (usually adopted) z reference 
direction is now the direction e3. This amounts to performing a rotation of the original variables 
Si to the new ones oj^ around a given axis and at a given angle depending on e3. Following the 
standard approach 



12l4l4j|. we use the Holstein-Primakov representation for the new variables <Tj 
. To rewrite the Green's functions in the local reference frame, we use a rotation matrix lZ(e y ,"d) 
for the rotation of an angle i? around the axis e y . So in the Hamiltonian ([T]) we replace the spin 
variable Sj by the new one <Xj (with ||<Xj|| = ||Sj|| = S) using 



Si = lZ(e y , — $)<Ti = cosher.; — sin?9 (cr, x e y ) + (1 — cos$) (<7j • e y )e y . (Al) 
For instance, 

S z = cos "daf - sin tfof 
and the Zeeman term (H y = 0) becomes 

Si H = [H x cos — H z sin <?] of + [H z cos •& + H x sin i?] af . ( A2) 

We also define the rotated field 

The new spin variables satisfy the same algebra as the original spin variables S$, i.e., 



26 tj af 



^■3 „M 

* > 3 



(A3) 
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Then, rewriting the Hamiltonian ([!]) in the new variables, we obtain the quadratic form 



H 



with the linear coefficients 



1 M N 

ij'=l|t,w=+,-,(3) i=l/i=+,-,(3) 



>)^, I? = {gv B ) Hi 



and the quadratic ones 



Qj+ = f sin 2 Sij = Q", Qtf = \ [Jij + K sin 2 = Q 



Q 33 = + 2K cos 2 Mi, Q?+ = Q± 3 = -K sin cos My = Q 3 ~ = Qr. 3 . 



(A4) 



(A5) 



(A6) 



These satisfy the symmetry relation = . 

Applying the RPA decoupling to a homogeneous ferromagnet, i.e. with (of) = (°" 3 )) we obtain 
the following (coupled) equations for the relevant GFs after Fourier transformations with respect 
to time and space 



where 



f uj-Ak 6 k 2A^ 

-B k u + A k -2A+ (K a ) 

K A+(JC a /2) -A k ()C a /2) u> 



( G+- \ ( 2 (a 3 ) \ 




/ 



ST 



Ak = L 3 + K a (a 3 ) (2 cos 2 - sin 2 0) + J (a 3 ) (1 - 7k ) , 
At = L± - /C (T (fj 3 }sin2i?, 
fik = £ CT (<r 3 ) sin 2 0. 
Jo is the k = component of the exchange coupling given by 



(A7) 



(A8) 



Jo = J (0) = Jij = Z J 

3 

with z being the coordination number. If the exchange is isotropic, we may write 
J(k) = J(-k) Y,' ' k V J U = J o x -E e " ik ' r " = J 07k- 



(A9) 



(A10) 
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For a bcc lattice we have (z = 8) the unit cell unit vectors 



and thereby 



Clkx QjJvy 0bl\>2, 

7k = cos cos — - cos — — , 

2 2 2 



(All) 



o being the lattice parameter. For long wavelength excitations we use cos/c a ~ 1 — ^k^, which 



yields 



1 - 7k - (akf 



Note that the EM for Q z - , that is the last equation in the system (JATJ) , provides the equilibrium 



configuration. Near equilibrium, the net magnetic moment m = i ^ i Cj does not change much, 
which means that dm/dt ~ 0. In quantum mechanics, this implies that the total magnetic moment 
along the equilibrium direction commutes with the Hamiltonian, or in other words, the projection 
of the total magnetic moment along the equilibrium direction is conserved, that is 



i dt\AfX' 



07 



0. 



(A12) 



On the other hand, on the same level of approximation as that used to obtain the system of EM 
(|A7j) . the commutator above reads 



\- r J2 a i' n ~ [(9Vb) (H z sm$ - H x cos$)+lC (a 3 ) sm2$] x -L ^ (ia?) , 



which, if set to zero according to (|A12p . leads to the equilibrium condition 



(h z sin d-h x cos •&) + K a (a 3 ) sin 2d = 0. 



(A13) 



(A14) 



Hence, the GF (ui) is eliminated from the system <\A7\i and thereby the latter simplifies into the 
following system of two coupled equations 

( 



2 a d 



0. 



(A15) 



Appendix B: Decoupling of exchange contributions within CSD 
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Following the procedure described in section UlI C[ before Eq. (j62[) . we obtain 
jf E 4 ({*- q WOO; C k (0)}) - mJ k ({4W;C k )) u 



TV 2 ^ q 



\ s p s q, 



'k-p-q 



p,q 



Similarly, for the second contribution we get 

Jf E 4 ({<(-)4- q W; s =k(0)}) w - mJj «4 (r); sZ k )) u 



— Tj' 



S ? S k-q 



(0 



q-pl T J; S -k 



p,q 



Now, using the two moment equations 

U-° oo^A k M = -<{{^,^},5i k }>. 

with the spectral density, see Ref. Q 

Ak M = «s+(r); S I k » w = i <{ S + (r) , sl k (0)})^ = 2*-jVm<J (w - a*) 
we integrate over u and obtain for the first contribution 



du 1 v . 

^AAE J ^i S k-ql^q 

q 



W;»: k (o)}) — 27Vm 2 j k - ^ ^ 



i \ s -q s q / 



and for the second 

oo 



-oo 1 1 

Then, subtracting the second contribution from the first yields 



k S k-q 



f v 

J q 



{*-,WM;-: k (o)}) 
{<(r)4- q (^);^k(o)} 



(Bl) 



(B2) 



2AAm 2 (J k -J^)-^^(j;-J k _ q ) 



\ S -q S q, 



On the other hand, using the zero-moment equation, we get 



oo 

f V 
J 2ttAT^ q 



{4- q W<(r); S : k (o)}) 
{«5(T)«t q (T);«: k (o)} 
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Consequently, we have the equation 



^ ^V"q "K-qy \ -q q / - ^ 

q q 



which leads to 

E K - J U) «*-q> - ^ 2 ™ 2 ( J o - JL)+\™ 2 E K " J U) <-iq-^ >4 E K - J k- q ) 

q q q 

One can can easily check that (s^ q s~) = (s! q Sq) and thereby one obtains 



E (4 " J U) «*-q> -EW- J k^) k- 2 A (q) - i (1 - m 2 ) < fl + fl = q > 
q q L 

This may also be recast into the form (|62l) which can be more easily compared to RPA. 



Appendix C: The critical temperature via different approaches. 

Within the QGF approach and using parameter (f> for exchange decoupling, the Curie tempera- 
ture can be calculated from Eq. (|43|) by setting (c 3 ) — at r = r c . This leads to [see Eq. ()33[) for 
notation] 

5(5 + 1)k + Q(^),t c ) 
Tc ~ 3 P(A) ' 

In the absence of anisotropy, which is negligible near T c , we obtain 



QGF = S(S + 1) Qexch(^,-r e QGF ) r QGF = 5(5+1) 

3 W c 3W 



^fl 1 + s) ( l -w 



(CI) 



In CGF (or in CSD where we obtain the same result), we similarly use the high-temperature 
asymptote (|60p and obtain 



; QU(<^ C CGF ) 



r CGF _ C 2 ^cxch . 

— O t T, 



CG¥\ Q 2 

CGF _ J 

c 



(C2) 



3W c 3W 

We remark in passing that this is also the result that one obtains within the spherical model, in 
the isotropic case j^l | , i.e. k = 0, and for a RPA decoupling cp = 0, which yields 
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On the other hand, from the MFT magnetization (|54p one obtains the Curie temperature (for 
H = and k = 0) 



MFT 



JqS 2 



(C4) 



Note that contrary to the MFT result (|C4p . the expression (|C2p for re, as obtained from the GF 
in the classical limit, or Eq. (|C3P from the isotropic spherical model, depends on the lattice and on 
the SW dispersion via the Watson integral W. Moreover, as mentioned earlier, we can relate MFT 
to SWT by assuming that all excitations are degenerate and by ignoring spin fluctuations. More 
precisely, this amounts to dropping the terms that are responsible for the propagation of the SWs 
(or magnons) through the lattice. This can be done by dropping the propagation function 7^ from 
all SWT expressions. Hence, the MFT result ()C4p can be obtained from the classical limit of the 
GF result (|C2p by formally setting 7k = in the lattice integral W (leading to W = 1) and taking 
= 0. 

In Table [I] we collect the values of tq estimated by the different approaches in the isotropic 
case. First, we remark that the values obtained within quantum-mechanical approaches are higher 
than the classical ones. Indeed, comparing for instance Eqs. (jCip and ()C2p we see that for small 
spin values the difference in r c , due to the contribution S in S (S + 1) = S 2 + S, is non negligible. 
This is no surprise because this decoupling scheme, unlike RPA, accounts for magnon-magnon 
interactions whose role becomes crucial in the vicinity of the critical temperature. Second, there 
is a perfect agreement between the two classical methods CGF and CSD. As discussed earlier, this 
shows that given that i) CGF renders the same results as CSD and ii) that CGF does not require 
any assumptions about the spectral function, it might be preferable to use the CGF method. 



Method 


QGF 

(a) (b) (c) 


CGF 

(a) (b) (c) 


CSD 

(a) (b) (c) 


MFT 


MC 


rc/S 2 (K) 


0.335 0.380 0.354 


0.240 0.262 0.262 


0.240 0.262 0.262 


0.333 


0.268 



Table I: Reduced Curie temperature tq/S 2 = k^Tc/ S 2 Jo for a bee lattice with Fe parameters and S = 5/2. 
(a) stands for the RPA or Copeland-Gersch, (b) for Callen, and (c) for Swenden's decoupling schemes used 
for the exchange contributions. 
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